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Abstract 


In  many  applications  of  tomography,  the  fundamental  quantities  of  interest  in  an  image  are  geometric  ones. 
In  these  instances,  pixel  based  signal  processing  and  reconstruction  is  at  best  inefficient,  and  at  worst,  non- 
robust  in  its  use  of  the  available  tomographic  data.  Classical  reconstruction  techniques  such  as  Filtered 
Back-Projection  tend  to  produce  spurious  features  when  data  is  sparse  and  noisy;  and  these  “ghosts”  further 
complicate  the  process  of  extracting  what  is  often  a  limited  number  of  rather  simple  geometric  features.  In 
this  paper  we  present  a  framework  that,  in  its  most  general  form,  is  a  statistically  optimal  technique  for 
the  extraction  of  specific  geometric  features  of  objects  directly  from  the  noisy  projection  data.  We  focus  on 
the  tomographic  reconstruction  of  binary  polygonal  objects  from  sparse  and  noisy  data.  In  our  setting,  the 
tomographic  reconstruction  problem  is  essentially  formulated  as  a  (finite  dimensional)  parameter  estimation 
problem.  In  particular,  the  vertices  of  binary  polygons  are  used  as  their  defining  parameters.  Under  the 
assumption  that  the  projection  data  are  corrupted  by  Gaussian  white  noise,  we  use  the  Maximum  Likelihood 
(ML)  criterion,  when  the  number  of  parameters  is  assumed  known,  and  the  Minimum  Description  Length 
(MDL)  criterion  for  reconstruction  when  the  number  of  parameters  is  not  known.  The  resulting  optimization 
problems  are  nonlinear  and  thus  are  plagued  by  numerous  extraneous  local  extrema,  making  their  solution  far 
from  trivial.  In  particular,  proper  initialization  of  any  iterative  technique  is  essential  for  good  performance. 
To  this  end,  we  provide  a  novel  method  to  construct  a  reliable  yet  simple  initial  guess  for  the  solution.  This 
procedure  is  based  on  the  estimated  moments  of  the  object,  which  may  be  conveniently  obtained  directly 
from  the  noisy  projection  data. 
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1  Introduction 


In  many  applications  of  tomography,  the  aim  is  to  extract  a  rather  small  set  of  geometrically-based  features 
from  a  given  set  of  projection  data  [1,  2,  3].  In  these  instances,  a  full  pixel-by-pixel  reconstruction  of  the 
object  is  a  rather  inefficient  and  non-robust  approach.  In  addition,  in  many  situations  of  practical  interest, 
a  full  set  of  data  with  high  signal-to-noise  ratio  (SNR)  is  often  difficult,  if  not  impossible,  to  obtain.  Such 
situations  arise  in  oceanography,  nuclear  medicine,  surveillance,  and  non-destructive  evaluation  when  due  to 
the  geometry  of  the  object  or  the  imaging  apparatus,  only  a  few  noisy  projections  are  available  [4,  5].  In 
these  cases,  the  classical  reconstruction  techniques  such  as  Filtered  Back  Projection  (FBP)  [5]  and  Algebraic 
Reconstruction  Techniques  (ART)  fail  to  produce  acceptable  reconstructions.  The  shortcomings  of  these 
classical  techniques  in  such  situations  can  be  attributed  to  two  main  sources.  First,  these  techniques  are 
invariably  aimed  at  reconstructing  every  pixel  value  of  the  underlying  object  with  little  regard  to  the  quality 
and  quantity  of  the  available  data.  To  put  it  differently,  there  is  no  explicit  or  implicit  mechanism  to  control 
greed  and  focus  information,  thus  preventing  one  from  attempting  to  extract  more  information  from  the 
data  than  it  actually  contains.  The  second  type  of  shortcoming  results  from  the  fact  that  if  we  assume 
that  the  projection  data  are  corrupted  by  Gaussian  white  noise,  the  process  of  reconstruction  will  have 
the  net  effect  of  “coloring”  this  noise.  This  effect  manifests  itself  in  the  object  domain  in  the  form  of 
spurious  features  which  will  complicate  the  detection  of  geometric  features.  This  observation  points  out  the 
importance  of  working  directly  with  the  projection  data  when  the  final  goal  is  the  extraction  of  geometric 
information.  In  our  effort  to  address  these  two  issues,  we  have  proposed  the  use  of  simple  geometric  priors 
in  the  form  of  finitely  parameterized  objects  (namely,  binary  polygonal).  The  assumption  that  the  object  to 
be  reconstructed  is  finitely  parameterized  allows  for  the  tomographic  reconstruction  problem  to  be  posed  as 
a  finite  (relatively  low-dimensional)  parameter  estimation  problem.  If  we  further  assume,  as  we  have  done 
in  the  latter  part  of  this  paper,  that  the  number  of  such  parameters  is  also  an  unknown,  we  can  formulate 
the  reconstruction  problem  as  a  Minimum  Description  Length  estimation  problem  which  provides  for  an 
automatic  (data-driven)  method  for  computing  the  optimal  parameterized  objects  with  the  “best”  number 
of  parameters,  given  the  data.  This  is,  in  essence,  an  information-theoretic  criterion  which  gives  us  a  direct 
way  to  estimate  as  many  parameters  as  the  information  content  of  the  data  allows  us  to,  and  thus  control 
the  greed  factor. 

Other  efforts  in  the  parametric/geometric  study  of  tomographic  reconstruction  have  been  carried  out  in 
the  past.  The  work  of  Rossi  and  Willsky  [6]  and  Prince  and  Willsky  [7,  8,  9]  has  served  as  the  starting 
point  for  this  research  effort.  In  the  work  of  Rossi,  the  object  was  represented  by  a  known  profile,  with  only 
three  geometric  parameters;  namely  size,  location,  eccentricity,  and  orientation.  These  parameters  were  then 
estimated  from  projection  data  using  the  Maximum  Likelihood  (ML)  formulation.  In  their  approach,  the 
number  of  unknown  parameters  was  fixed  and  the  main  focus  of  their  work  was  on  performance  analysis. 
Prince,  on  the  other  hand,  used  a  priori  information  such  as  prior  probabilities  on  sinograms  and  consistency 
conditions  to  compute  Maximum  A  Posteriori  (MAP)  estimates  of  the  sinogram  and  then  used  FBP  to 
obtain  a  reconstruction.  He  made  use  of  prior  assumptions  about  shape,  such  as  convexity,  to  reconstruct 
convex  objects  from  support  samples  which  were  extracted  from  noisy  projections  through  optimal  filtering 
techniques.  The  approach  of  Prince  provided  an  explicit  method  for  integrating  geometric  information  into 
the  reconstruction  process  but  was  in  essence  still  a  pixel-by-pixel  reconstruction.  Extending  these  ideas,  Lele, 
Kulkarni,  and  Willsky  [10,  11]  made  use  of  only  support  information  to  produce  polygonal  reconstructions. 
Hanson  [12]  studied  the  reconstruction  of  axially  symmetric  objects  from  a  single  projection.  Karl  [13]  also 
has  studied  the  reconstruction  of  3-D  objects  from  two-dimensional  silhouette  projections. 

The  geometric  modeling  approach  of  Rossi  and  Willsky  was  expanded  upon  to  include  a  more  general 
set  of  objects  by  Bresler,  Macovski  and  Fessler  [14,  15,  16,  17].  In  these  papers,  the  authors  chose  sequences 
of  3-D  cylinders  and  ellipsoids  parameterized  by  stochastic  dynamic  models  based  on  their  radius,  position, 
and  orientation  to  model  and  reconstruct  objects  in  two  and  three  dimensions  from  projections. 

Recently,  Thirion  [2]  has  introduced  a  technique  to  extract  boundaries  of  objects  from  raw  tomographic 
data  through  edge  detection  in  the  sinogram.  Other  work  in  geometric  reconstruction  by  Chang  [18]  and  more 
recently  Kuba,  Volcic,  Gardner  and  Fishburn,  [19,  20,  21,  22]  has  been  concerned  with  the  reconstruction  of 
binary  objects  from  only  two  noise-free  projections. 

Our  approach  provides  a  statistically  optimal  ML  formulation  for  the  direct  recovery  of  vertices  of  binary 
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Figure  1:  The  Radon  transform 


polygons  from  the  projection  data  in  the  presence  of  noise.  We  also  provide  an  automatic  mechanism  for 
identifying  the  statistically  optimal  number  of  vertices,  from  a  given  data  set.  The  statistically  optimal  ML 
formulation  leads  to  an  optimization  problem  that  is  nonlinear  and  filled  with  local  extrema.  An  appropriate 
initial  guess  is  thus  essential  for  its  iterative  solution.  An  important  contribution  of  this  paper  is  that  we 
thus  provide  a  simple  procedure  to  generate  an  appropriate  initial  guess  based  on  moment  estimates  of  the 
object  computed  from  the  original  projection  data. 

The  organization  of  this  paper  is  as  follows.  In  Section  2  we  introduce  the  basic  definitions  and  as¬ 
sumptions  and  pose  the  general  problem  which  we  intend  to  solve.  We  also  discuss  the  particular  statistical 
formulations  of  the  reconstruction  problem  which  we  use.  In  particular,  in  Section  2.3  we  discuss  our  novel 
technique  for  computing  a  good  initial  guess  for  the  nonlinear  optimization  problems  that  result  from  our 
formulations.  Section  3  contains  basic  performance  results  and  robustness  studies  for  various  scenarios. 
Section  4  contains  our  conclusions. 


2  The  Reconstruction  Problem 

The  Radon-Transform  [5,  23]  of  a  function  f(x,  y)  defined  over  a  compact  domain  of  the  plane  O  is  given  by 

g(t,  6)  =  J f(x,  y)  S(t-u-  [x,  y]T )  dxdy.  (1) 

For  every  fixed  t  and  9,  g(t,  9)  is  the  line-integral  of  /  over  O  in  the  direction  to  =  [cos(0),  sin(0)]T,  where 
S(t  —  [cos(0),  sin(0)]  •  [x,  y)T)  is  a  delta  function  on  a  line  at  angle  9  +  7r/2  and  distance  t  from  the  origin. 
See  Figure  1. 

Here  we  assume  that  the  function  /  is  the  indicator  function  for  some  simply-connected  binary  polygon 
and  hence  a  finite  set  of  parameters  uniquely  specify  the  function  /.  The  estimation  of  the  parameters  that 
uniquely  specify  the  function  /  is  the  concern  of  this  paper.  Let  us  stack  the  polygon  vertices  that  uniquely 
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define  /  in  a  vector  V.  We  will  assume  throughout  that  the  data  available  to  us  are  discrete  samples  of 
g  which  are  corrupted  by  Gaussian  white  noise  of  known  intensity.  In  particular,  our  observations  will  be 
given  by 

Yij  =g(ti,ej,V*)  +  w(tuej),  (2) 

for  1  <  i  <  m,  1  <  j  <  n  where  V*  is  the  true  object  we  wish  to  reconstruct.  The  variables  w(ti,0j )  are 
assumed  to  be  independent,  identically  distributed  (i.i.d.)  Gaussian  random  variables  with  variance  er2.  We 
will  denote  by  Y  the  vector  of  all  such  observations. 

2.1  Maximum  Likelihood  Approach 

In  our  approach,  the  original  data  in  (2)  are  used  to  directly  estimate  the  vertices  V  of  a  polygon  in  a 
statistically  optimal  way.  The  dimension  of  the  parameter  vector  V  (i.e.  the  number  of  sides)  is  determined 
by  the  level  of  detail  that  one  can  extract  from  the  sparse  and  noisy  data.  For  clarity,  we  first  consider  the 
case  where  a  fixed  and  known  number  of  vertices  is  assumed.  In  this  case,  the  Maximum  Likelihood  (ML) 
[24]  estimate,  Vmi,  of  the  parameter  vector  V  is  given  by  that  value  of  V  which  makes  the  observed  data 
most  likely.  In  particular,  using  the  monotonicity  of  the  logarithm: 

Vmi  =  argmax  log  [P(Y\V)] ,  (3) 

where  P(Y\V)  denotes  the  conditional  probability  density  of  the  observed  data  set  Y  given  the  parameter 
vector  V.  It  is  well-known  that  given  the  assumption  that  the  data  is  corrupted  by  i.i.d.  Gaussian  random 
noise,  the  solution  to  the  above  ML-estimation  problem  is  precisely  equivalent  to  the  following  Nonlinear 
Least  Squares  Error  (NLSE)  formulation 

Vmi  =  argmin^  || Y{J  -  g(ti,0j,V) ||2.  (4) 

i,j 

The  formulation  (4)  shows  that,  in  contrast  to  the  linear  formulation  of  classical  reconstruction  algorithms, 
the  ML  tomographic  reconstruction  approach,  while  yielding  an  optimal  reconstruction  framework,  generally 
results  in  a  highly  nonlinear  minimization  problem.  It  is  the  nature  of  the  dependence  of  g  on  the  parameter 
vector  V  that  makes  the  problem  nonlinear. 

Finally,  note  that  if  additional  explicit  geometric  information  is  available  in  terms  of  a  prior  probabilistic 
description  of  the  object  vector  V,  then  a  Maximum- A-Posteriori  estimate  of  V  may  be  computed  as  follows: 

Vmap  =  argmax  log  [P(V\Y)}  (5) 

In  this  work  we  concentrate  on  the  ML  problem  given  in  (3)  and  its  extensions,  though  application  of  our 
results  to  the  MAP  formulation  is  straightforward. 

2.2  Minimum  Description  Length 

In  the  previous  ML  discussion  we  assumed  we  had  prior  knowledge  of  the  number  of  parameters  describing 
the  underlying  object.  Without  this  knowledge,  we  can  consider  the  Minimum  Description  Length  (MDL) 
principal  [25].  In  this  approach,  the  cost  function  is  formulated  such  that  the  global  minimum  of  the  cost 
corresponds  to  a  model  of  least  order  that  explains  the  data  best.  The  MDL  approach  in  essence  extends  the 
Maximum  Likelihood  principal  by  including  a  term  in  the  optimization  criterion  that  measures  the  model 
complexity.  In  the  present  context,  the  model  complexity  refers  to  the  number  of  vertices  used  to  capture 
the  object  in  question.  Whereas  the  ML  approach  maximizes  the  log  likelihood  function  given  in  (3),  the 
MDL  criterion  maximizes  a  modified  log  likelihood  function,  as  follows: 

Vmdl  =  argmax  jlog  [P(Y\ V)}  -  y  log (d)  j  ,  (6) 

where  d  =  mn  is  the  number  of  samples  of  g(t,  6)  and  N  refers  to  the  number  of  parameters  defining  the 
reconstruction.  Roughly  speaking,  the  MDL  cost  is  proportional  to  the  number  of  bits  required  to  model 
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Figure  2:  A  projection  of  a  binary,  polygonal  object 


the  observed  data  set  with  a  model  of  order  N:  hence  the  term  Minimum  Description  Length.  Under  our 
assumed  observation  model  (2)  the  MDL  criterion  (6)  yields  the  following  nonlinear  optimization  problem 
for  the  optimal  parameter  vector  Vmdl • 

Vmdl  =  ar  g  min  min  { <7 -  2  ^  ||  Yitj  -  g(U,  Oj,V)  ||2  +  TVlog(d)},  (7) 

i,j 

where  the  optimization  is  now  performed  over  both  V  and  the  number  of  parameters  N.  Note  that  the 
solution  of  the  inner  minimization  in  (7)  essentially  requires  solution  of  the  original  ML  problem  (3)  or  (4) 
for  a  sequence  of  values  of  N. 

Unless  otherwise  stated,  we  assume  from  here  on  that  the  matrix  V  contains  the  vertices  of  an  iV-sided 
binary,  polygonal  region  as  follows: 

V  =  [  UL  t  V2  I  •  •  •  I  VN  ]  ,  (8) 

where  Vj  =  [xi,yi]T  denote  the  Cartesian  coordinates  of  the  ith  vertex  of  the  polygonal  region  arranged  in 
the  counter-clockwise  direction  (See  Figure  2).  Note  that  we  use  a  matrix  of  parameters  rather  than  a  vector 
in  what  follows  for  notational  convenience  in  the  algorithms  to  follow,  though  this  is  not  essential. 

2.3  Algorithmic  Aspects  —  Computing  A  Good  Initial  Guess 

Given  the  highly  nonlinear  nature  of  the  dependence  of  the  cost  function  in  (4)  on  the  parameters  in  V ,  it 
appears  evident  that  given  a  poor  initial  condition,  typical  numerical  optimization  algorithms  may  converge 
to  local  minima  of  the  cost  function.  Indeed,  this  issue  is  a  major  obstacle  to  the  use  of  a  statistically  optimal, 
though  nonlinear,  approach  such  as  given  in  (3)  or  (6).  In  this  section  we  describe  a  method  for  using  the 
projection  data  to  directly  compute  an  initial  guess  that  is  sufficiently  close  to  the  true  global  minimum  as 
to,  on  average,  result  in  convergence  to  it,  or  to  a  local  minimum  nearby.  We  do  this  by  estimating  the 
moments  of  the  object  directly  from  the  projection  data  and  then  using  (some  of)  these  moments  to  compute 
an  initial  guess. 

In  considering  the  use  of  moments  as  the  basis  for  an  initialization  algorithm,  one  is  faced  with  two 
important  issues.  The  first  is  that  although  estimating  the  moments  of  a  function  from  its  projections  is 
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a  relatively  easy  task,  as  we  have  shown  in  [26,  27],  the  reconstruction  of  a  function  from  a  finite  number 
of  moments  is  in  general  a  highly  ill-posed  problem  even  when  these  moments  are  exactly  known  [28]. 
Furthermore,  in  our  framework  the  moments  are  estimated  from  noisy  data,  and  hence  are  themselves  noisy. 
In  fact,  as  higher  and  higher  order  moments  are  estimated,  the  error  in  the  estimates  of  these  moments 
becomes  larger.  Our  approach  avoids  these  moment  related  difficulties  by  using  the  moments  only  to  guide 
an  initial  coarse  estimate  of  the  object  parameters  for  subsequent  use  in  solution  of  the  nonlinear  ML  or  MDL 
problems.  This  initial  estimate,  in  turn,  itself  mitigates  the  difficulties  associated  with  the  nonlinearities  of 
the  optimal  statistical  approaches.  In  particular,  the  amount  of  computation  involved  in  arriving  at  an 
initial  guess  using  our  moment-based  method  is  far  smaller  than  the  amount  of  computation  (number  of 
iterations)  required  to  converge  to  an  answer  given  a  poor  initial  guess,  especially  since  a  poor  initial  guess 
may  converge  to  a  local  minimum  far  from  the  basin  of  the  global  minimum.  Further,  the  parameterization 
of  the  objects  serves  to  regularize  and  robustify  the  moment  inversion  process  [28,  29,  30,  31]. 

Our  method  of  using  moments  to  generate  an  initial  guess  is  based  on  the  following  set  of  observations. 
First,  let  y pq ,  0  <  p,q  denote  the  moment  of  f(x,  y)  of  order  p  +  q  as  given  by: 


Ppq 


xpyqf(x,y)  dxdy 


(9) 


In  particular,  note  that  the  moments  up  to  order  2  have  the  following  physical  relationships.  The  zeroth 
order  moment  yoo  is  the  area  of  the  object,  the  first  order  moments  yoi  and  ^io  are  the  coordinates  of  the 
center  of  mass  of  the  object  scaled  by  the  area,  and  the  second  order  moments  yo2,  Mn>  M20  are  used  to  form 
the  entries  of  the  inertia  matrix  of  the  object.  Thus  these  moments  contain  basic  geometric  information 
about  object  size,  location,  and  elongation  and  orientation  that,  if  available,  could  be  used  to  guide  our 
initialization  of  the  nonlinear  optimization  problems  (4)  or  (7).  Our  first  aim  then  is  to  estimate  them 
directly  from  the  noisy  projection  data.  To  that  end,  it  is  easy  to  show  that  [23]: 

g(t,9)tkdt  =  j  J  f(x,y)[xcos(6)+ysin(6)]kdxdy.  (10) 


By  expanding  the  integrand  on  the  right  hand  side  of  (10),  it  becomes  apparent  that  the  moments  of  the 
projections  are  linearly  related  to  the  moments  ypq  of  the  object.  In  particular,  specializing  (10)  to  k  =  0, 1,  2, 
and  noting  that  f(x,y)  is  an  indicator  function  when  the  objects  in  question  are  binary,  we  arrive  at  the 
following  relationships  between  ypq,  0  <  p  +  q  <  2  and  the  projection  g(t,  9)  of  f(x,  y)  at  each  angle  9: 


Moo 


cos  (9)  sin(0) 


cos 2(9)  2sin(0)  cos(M)  sin2(0) 


M10 

M01 

M20 

Mil 

M02 


J  g(t,9)dt  =  H^(9)  (11) 

I  g(t,9)tdt  =  nw(9)  (12) 

I  g(t,9)t2  dt  =  H(2)(9)  (13) 


Thus  if  we  have  projections  at  three  or  more  distinct,  known  angles  we  can  estimate  the  moments  of  up  to 
order  2  of  the  object  we  wish  to  reconstruct.  The  computation  of  these  moments  is  a  linear  calculation, 
making  their  estimation  from  projections  straightforward  (see  [26]).  Since,  in  general,  many  more  than  three 
projections  are  available,  the  estimation  of  these  moments  determining  the  area,  center,  and  inertia  axes  of 
the  object  is  overdetermined.  The  result  is  a  robustness  to  noise  and  data  sparsity  through  a  reduction  in 
the  noise  variance  of  their  estimated  values.  In  particular,  we  can  stack  the  moments  (9j)  obtained  from 
the  projections  at  each  angle  9:j  to  arrive  at  the  following  overall  equations  for  the  ypq  up  to  order  2: 


'  1  ' 

'  n(0\9i)  ' 

Moo  — 

1 

_  HW(9n)  . 

(14) 
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COS  (0i ) 

sin  (0i ) 

'  W(1>(0i)  ' 

Mio 

= 

_  cos (0„) 

sin(0„) 

M  oi 

.  «(')(«„)  . 

cos2  (6*i ) 

2sin(0i) 

cos(0i) 

sin2  (0i ) 

M20 

'  (0i )  ' 

Mil 

= 

cos2(0n) 

2sin(0„) 

cos(0n) 

sin2(0„)  _ 

M02 

.  H(2\0n)  _ 

(15) 


(16) 


Using  these  equations  we  can  easily  calculate  the  least  squares  (LS)  estimates  of  the  moments  of  the 
object  [ipq  for  0  <  p  +  q  <  2  given  noisy  observations  of  the  moments  of  the  projections.  In  par¬ 
ticular,  this  is  done  by  gathering  the  above  sets  of  equations  into  a  single  linear  equation  of  the  form 
h  =  A/i  +  e,  where  h  is  the  vector  of  noisy  computed  moments  of  the  projections  appearing  above, 

and  /z  =  [/zoo,  Mien  Moi>  M20>  Mil,  M02]2",  while  e  denotes  a  zero-mean  vector  of  Gaussian  noise  with  correspond¬ 
ing  covariance  matrix  R  which  captures  the  noise  in  our  observations  of  TL(k'>  (0,).  This  noise  model,  of  course, 
follows  directly  from  that  of  (2).  The  least  squares  estimate  of  the  vector  /z  is  then  given  by 


Ji  =  (ATR-1A)-1ATR~1h. 


(17) 


Note  that  if  this  least  squares  estimate  is  consistent1,  it  will  coincide  with  the  optimum  Maximum  Likelihood 
(ML)  estimate  of  the  moments  of  the  underlying  polygon.  The  estimates  of  the  low  order  moments,  such 
as  those  we  use  here,  are  robust  to  noise  in  the  data  and  hence  in  the  great  majority  of  the  cases,  the  LS 
estimate  does  indeed  coincide  with  the  ML  estimate.  When  this  fails,  it  is  typically  due  to  the  inconsistency 
of  the  second  order  moment  estimates.  In  such  cases,  we  refrain  from  using  the  second  order  estimates 
and  build  our  initial  guess  using  only  the  area  and  center  of  mass  estimates  as  we  describe  next.  (The 
general  framework  for  the  optimal  estimation  of  the  moments  of  any  order  of  a  function  /(a:,  y)  from  noisy 
measurements  of  its  Radon  transform  is  developed  in  [26,  27].)  We  shall  henceforth  denote  the  LS  moment 
estimates  by  /zP9. 

Now  that  we  have  estimates  of  the  moments  of  up  to  order  2  of  the  object,  and  thus  estimates  of  its 
basic  geometric  structure,  we  need  to  convert  this  information  into  a  suitable  object  for  use  in  initializing  the 
nonlinear  optimization  problem  (4)  or  (7).  The  initial  guess  algorithm  outlined  next  uses  these  least  squares 
estimates  of  the  low-order  moments  jj.pq ,  obtained  from  the  noisy  projection  data,  to  produce  a  polygon 
which  has  moments  up  to  order  2  which  are  close  to  (or  in  some  cases  equal  to)  those  which  were  estimated 
from  the  projection  data.  The  resulting  polygon,  which  will  be  used  as  our  initialization,  should  thus  have 
the  same  basic  geometric  structure  as  the  underlying  object. 

Recall  that  in  this  process  of  generating  an  initial  object  from  the  moment  data  we  want  to  avoid  the 
difficulties  usually  associated  with  obtaining  an  object  from  a  set  of  its  moments  [29,  31,  30,  28].  For  this 
reason,  the  initial  polygon  we  will  use  is  simply  obtained  as  the  affine  transformation  of  a  reference  object 
Vi.ef (TV),  which  is  a  centered  regular  N- gon  of  unit  area.  For  a  given  choice  of  number  of  sides  N,  the 
reference  object  we  use  is  given  by: 


Kef  (AO  = 


1 


cos(0)  cos(^) 
sin(0)  sin(^) 


cos (2=^)  ■ 
sin(M^il)  _ 


(18) 


The  affine  transformation  of  this  reference  object,  which  will  be  generated  from  the  estimated  moment  set, 
consists  of  a  uniform  scaling,  a  stretching  along  the  coordinate  axes,  a  rotation,  and  finally  a  translation, 
and  simply  corresponds  to  the  following  transformation  of  the  underlying  spatial  coordinates  of  the  reference 
object 


'  a/  ' 

=  L 

X 

L  y  \ 

.  y . 

+  C. 


(19) 


1that  is,  it  satisfies  the  necessary  and  sufficient  conditions  to  be  the  moment  vector  of  some  binary  polygon,  e.g.  give  a 
positive  area  estimate  and  a  positive  definite  inertia  matrix  estimate.  We  have  dealt  with  the  inconsistent  moment  case  directly 
in  [26]. 
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In  particular,  given  the  form  of  Vref(N)  in  (18),  this  yields  the  following  equation  for  the  family  of  possible 
initial  objects  Vjnit : 

Unit  =LVref(N)+  [  C  |  •••  |  C  ]  .  (20) 

The  set  of  all  such  affine  transformations  of  Vlef(N)  we  term  the  affinely  regular  IV-gons  [32].  In  the 
absence  of  noise,  the  initial  guess  algorithm  we  detail  will  exactly  match  the  given  estimated  moments  if  the 
underlying  object  itself  happens  to  be  affinely  regular.  If  the  underlying  object  is  not  affinely  regular  itself, 
the  algorithm  will  not  necessarily  produce  an  IV-gon  exactly  matching  its  moments,  even  in  the  absence  of 
noise,  though  as  we  will  show,  it  will  usually  be  close.  Of  course,  in  the  presence  of  noise  the  estimated 
moments  themselves  are  not  exact  and  thus,  while  we  would  hope  to  get  close,  our  resulting  initial  V-gon 
will  never  exactly  match  the  true  moments  of  the  underlying  object  anyway. 

Given  that  we  will  restrict  ourselves  to  initial  objects  of  the  form  (20),  let  us  consider  how  we  might  choose 
the  parameters  of  the  transformation  L  and  C  to  match  a  given  estimated  moment  set  0  <  p  +  q  <  2. 
Using  (20)  and  (18)  to  calculate  the  moments  of  Unit  up  to  order  2  we  obtain  the  following  relationships: 


Moo  (Unit) 
MlO  (Unit) 
M01  (Unit) 

M20  (Unit )  Mil  (Unit) 
Mil  (Unit)  M02  (Unit) 

where  pvq  (Unit)  is  the  pq- th  moment  of  Unit  and  fcjv 
Moi(Unit)  with  their  estimated  counterparts,  the  fi: 


det(L) 

(21) 

det(L)  C 

(22) 

det(L)  (kNLLT  +  CCT) 

(23) 

l/(4iVtan(7r/JV)).  Thus  to  match  Moo  (Unit),  Mio(Unit), 
two  conditions  require: 


|det(T)| 

C 


Moo 

1 

Moo 


M10 

M01 


(24) 

(25) 


The  first  condition  simply  corresponds  to  a  scaling  of  Uef  (IV)  so  that  its  area  matches  the  estimated  one.  The 
second  condition  shows  that  the  affine  term  C  in  the  transformation  (19)  should  correspond  to  a  translation 
of  Vief(N)  to  the  estimated  center  of  mass  of  the  object.  These  two  conditions  assure  that  we  match  the 
estimated  area  and  center  of  mass  location. 

Now,  after  some  manipulation,  (23)  implies  that  to  match  the  estimated  second  order  moments  p,pq, 
p  +  q  =  2,  we  must  have: 

LLt  =  -i- 1  (26) 

UvMoo 

where  X  is  the  matrix  of  estimated  central  moments  defined  by: 


X  = 

M20 

Mu 

1 

M10 

Mil 

M02 

Poo 

MO! 

M01 


(27) 


In  particular,  this  condition  implies  another  constraint  on  clet(T)  independent  of  (24),  which  we  will  not,  in 
general,  be  able  to  satisfy.  Specifically,  a  necessary  condition  for  finding  an  L  satisfying  both  (26)  and  (24) 
is  that: 

det(T)  =  fcjvMoo-  (28) 

Actually,  as  shown  in  Appendix  A,  the  condition  (28)  is  also  sufficient.  Clearly,  this  condition  will  not, 
in  general,  be  satisfied  and  we  will  be  unable  to  exactly  match  the  estimated  second  moments.  In  fact, 
the  objects  that  do  meet  this  constraint,  and  thus  whose  moments  we  can  exactly  match,  are  precisely  the 
elements  of  the  set  of  affinely  regular  N- gons.  Geometrically,  this  situation  reflects  the  limitation  of  our 
restricted  initial  guess  object  class  (20),  i.e.  the  set  of  affinely  regular  N- gons.  Within  this  class,  for  a  given 
object  area  we  are  constrained  as  to  the  “size”  of  the  corresponding  inertia  matrix  we  may  have,  where 
inertia  size  is  measured  as  the  product  of  the  principal  central  inertia  moments  (eigenvalues  of  the  central 
inertia  matrix).  For  example,  while  our  initial  guess  objects  will  always  be  convex  polygons,  for  a  given  area 
we  can  always  obtain  nonconvex  objects  of  greater  inertia  by  “moving  area  outward,”  as  in  a  dumbell. 
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The  condition  (26)  can  also  be  viewed  as  implying  a  different  scaling  on  L  needed  to  obtain  a  perfect 
match  to  the  inertia  condition.  In  general,  we  thus  have  a  choice  of  picking  this  scaling  of  L  to  satisfy  either 
the  inertia  condition  (26)  or  the  area  condition  (24).  Since  the  area  condition  (24)  is  both  a  more  robustly 
estimated  and  a  more  fundamental  geometric  quantity,  we  choose  to  enforce  this  condition  in  the  algorithm 
to  follow.  We  then  choose  L  so  that  the  resulting  central  inertia  matrix  of  Vjnit  has  the  same  principal  axes 
directions  and  has  its  principal  inertias  in  the  same  ratio  as  those  estimated  from  the  data  as  found  in  X.  We 
accomplish  these  goals  by  using  a  square  root  of  the  matrix  X  normalized  to  have  unit  determinant  for  the 
form  of  L ,  then  scaling  the  result  to  match  the  determinant  condition  (24).  Thus  we  sacrifice  overall  scaling 
of  the  inertia  matrix  in  favor  of  matching  the  estimated  area.  Collecting  the  above  steps  and  reasoning,  the 
overall  algorithm  is  given  by  the  following: 

Algorithm  1  (Initial  Guess  Algorithm) 

1.  Compute  the  least  squares  estimates  of  the  moments  up  to  order  2  (fioo,  Mio.  Moi.  M20<  M11,  and  P02) 
from  the  raw  projection  data  using  (lf)-(l  7). 

2.  Construct  an  N -sided  regular  polygon  centered  at  the  origin  with  vertices  chosen  as  the  scaled  roots  of 
unity  in  counter-clockwise  order  so  that  they  lie  on  the  circle  of  radius  1/ yj sin(^).  This  polygon 
has  unit  area  and  is  defined  in  Equation  18. 

3.  Compute  the  translation  C,  obtained  as  the  estimated  object  center  of  mass: 


1  M 10 
Moo  .  M01 


(29) 


4-  Form  the  estimated  central  inertia  matrix  X  from  the  estimated  moments  according  to  (27) 

5.  IfX  is  not  positive  definite,  set  L  =  a/moo  I2  and  goto  step  8.  Otherwise  proceed  to  step  6.  (I2  is  the 
2x2  identity  matrix.) 

6.  Perform  an  eigendecomposition  of  the  normalized  matrix  X  as  follows: 

=  Ul«  1/A  ]  ^  (30) 

where  we  have  assumed  that  the  eigenvalues  are  arranged  in  descending  order  and  that  the  eigenvectors 
are  normalized  to  unit  length  so  that  det(CZ)  =  ±1.  Note  that  the  eigenvalues  are  reciprocals  of  each 
other  since  we  have  scaled  the  left,  hand  side  so  that  its  determinant  is  1. 

7.  Form  the  linear  transformation  L  as  a  scaled  square  root  of  X  as  follows: 

L  =  ^ 17  [  f  l/VX  ]  <31) 

Note  that  det(T)  =  moo  as  desired.  Depending  on  whether  U  is  a  pure  rotation  or  a  rotation  followed 
by  a  reflection,  it  will  have  determinant  +1  or  — 1,  respectively.  Also  note  that  we  use  this  (Schur 
decomposition-based)  square  root  formulation  over  the  more  numerically  attractive  Cholesky  factor¬ 
ization  since  our  choice  leads  to  matrices  that  are  directly  interpretable  as  rotations,  reflections,  and 
scalings. 

8.  The  initial  guess  V)nit  is  now  obtained  by  applying  the  scaling,  stretching,  and  rotation  transformation 
L  and  the  translation  C  to  the  reference  object  Vre{ (N)  via  the  coordinate  transformation  [x' ,  y')T  = 
L[x,  y)T  +  C.  Because  of  the  form  ofVie{(N )  this  operation  yields: 

Mnit  =  LViet(N)  +  [  C  |  •••  |  C  ]  .  (32) 
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Note  that  the  eigenvalue  A  of  the  unit  determinant  matrix  calculated  in  step  7  gives  the  eccentricity  of  the 
underlying  object  while  the  corresponding  eigenvectors  give  its  orientation.  Also  note  that  in  the  presence 
of  noise  the  estimated  central  inertia  matrix  for  the  object  X  may  not  be  strictly  positive  definite  and  hence 
may  not  correspond  to  the  inertia  matrix  of  any  object  at  all.  In  such  instances,  the  algorithm  refrains  from 
the  use  of  these  moments  of  order  2  and  computes  an  initial  guess  based  only  on  the  estimated  area  and 
center  of  mass. 

If  the  matrix  L ,  computed  above,  is  replaced  by  L'  =  LT  for  any  orthogonal  T,  the  resulting  quantity 
L'L't  satisfies 

LLt  =  L’L’t.  (33) 

Hence,  although  the  initial  guess  generated  by  the  above  algorithm  is  unambiguous  and  unique  in  the  sense 
that  the  square  root  of  X  obtained  by  the  algorithm  is  unique,  an  infinity  of  other  initial  guesses  having  the 
same  moments  up  to  order  2  may  be  generated  by  replacing  L  by  LT  and  allowing  T  to  range  over  the  set 
of  all  2  x  2  orthogonal  transformations.  A  precise  characterization  of  this  set  of  all  affinely  regular  IV-gons 
with  the  same  moments  up  to  order  2  is,  in  fact,  given  in  the  following  result,  which  we  prove  in  Appendix 
A. 


Result  1  Consider  the  set  of  all  N-gons  with  moments  up  to  order  2:  /zoo.  Mio>  Moi>  M20.  /zii,  M02.  such 
that  the  resulting  inertia  matrix 


M  20  M11 

M11  M  02 


(34) 


satisfies  det(X)  =  fc^Moo-  This  set  coincides  with  the  set  of  N-gons  with  vertices  on  an  ellipse  £ o  and  sides 
tangent,  at  their  midpoints,  to  a  homofocal  ellipse  £1,  where  these  ellipses  are  given  by 


£o  =  {z\(z-C)tEo1(z-C)  =  1} 
£I  =  {z\(z-C)TEj\z-C)  =  l}. 


with 


c  = 

1 

Moo 

M10 

M01 

Eo  = 

4 

/Zoo  cos2  (71-/ A) 

Ei  = 

— X 

Moo 

(35) 

(36) 

(37) 

(38) 

(39) 


This  result  states  that  the  class  of  all  affinely  regular  A-gons  with  a  given  set  of  moments  up  to  order  2  is 
given  by  the  class  of  A-gons  whose  vertices  are  on  a  fixed  ellipse  and  whose  side  are  tangent  to  a  second 
ellipse  which  is  homofocal  with  the  first.  (See  Figure  3  for  an  example.)  The  ellipses  are  uniquely  determined 
by  the  value  of  the  given  moments.  The  above  result  draws  attention  to  a  question  of  a  more  general  nature. 
Namely,  “How  many  moments  (and  hence  projections)  uniquely  specify  a  simply-connected  Ar-gon?”  We 
answer  this  question  in  [33]  and  [26]  where  in  the  former  we  have  shown  that  moments  up  through  order 
2A  —  3  suffice  to  uniquely  specify  the  vertices  of  any  simply-connected  A-gon,  while  in  the  latter  we  have 
shown  that  m+  1  projections  are  necessary  and  sufficient  to  uniquely  specify  the  moments  up  through  order 
to  of  any  object  from  its  projections.  Hence,  together  these  results  show  that  with  2N  —  2  projections,  we 
can  uniquely  specify  the  vertices  of  any  A-sided,  simply-connected  polygon.  In  contrast  to  [33,  26]  where 
image  reconstruction  from  projections  is  based  directly  on  the  estimated  moments,  the  initial  guess  algorithm 
presented  here  aims  only  to  give  a  rough  initial  guess  which  is  hopefully  within  a  reasonable  neighborhood 
of  the  global  optimum  of  the  ML  cost. 

In  order  to  simplify  the  Initial  Guess  Algorithm,  we  do  not  search  further  over  the  family  Given  by  Result 
1.  We  simply  use  the  output  of  the  Initial  Guess  Algorithm  described  above  as  the  starting  guess  for  our 
nonlinear  optimization  routines. 


13 


Figure  3:  Illustration  of  Result  1 


3  Experimental  Results 

In  this  section  we  present  some  performance  studies  of  our  proposed  algorithm  with  simple  polygonal  objects 
as  prototypical  reconstructions.  One  may  expect  that  our  algorithms  work  best  when  the  underlying  object 
(that  from  which  the  data  is  generated)  is  itself  a  simple  binary  polygonal  shape.  While  this  is  true,  we  will 
also  show  that  our  algorithms  perform  well  even  when  the  underlying  objects  are  complex,  non-convex,  and 
non-polygonal  shapes. 

First  we  demonstrate  reconstructions  based  on  the  ML  criterion.  In  these  reconstructions  we  use  the 
parameters  of  the  true  polygon  as  the  initial  guess.  Typical  reconstructions  are  shown  along  with  average 
performance  studies  for  a  variety  of  noise  and  data  sampling  scenarios.  In  particular,  we  show  that  given 
a  good  initial  guess,  the  performance  of  our  algorithms  is  quite  robust  to  noise  and  sparsity  of  the  data, 
significantly  more  so  than  classical  reconstruction  techniques.  To  demonstrate  this  point,  reconstructions 
using  our  techniques  and  the  classical  FBP  are  provided. 

Next  we  demonstrate  how  the  MDL  criterion  may  be  used  to  optimally  estimate  the  number  of  parameters 
(sides)  N  directly  from  the  data.  We  solve  these  MDL  problems  by  solving  the  ML  problem  for  a  sequence  of 
values  of  N.  To  initialize  each  of  these  ML  problems,  the  Initial  Guess  Algorithm  was  used.  The  robustness 
of  the  MDL  approach  and  its  ability  to  capture  the  shape  information  in  noisy  data  when  the  underlying 
object  is  not  polygonal  is  also  shown  through  polygonal  reconstruction  of  more  complicated  shapes.  Two 
remarks  are  in  order  regarding  the  results  presented  in  Section  3.1.  The  first  is  that  the  results  in  this 
sections,  begin  by  initializing  the  algorithms  with  the  tru  underlying  polygon.  This  correspond  to  local 
exploration  of  the  ML  and  MDL  cost  functions.  Hence,  any  interpretation  of  these  performance  results  (e.g. 
robustness)  is  conditioned  on  whether  or  not  the  initial  guess  algorithm  will  provide  a  starting  guess  that 
is  within  the  local  basin  of  the  global  optimum  (which  as  we  show  in  Section  3.1  apparently  includes  the 
underlying  polygon) .  In  the  subsequent  sections  we  show  that  the  Initial  Guess  Algorithm  does,  on  average 
place  the  starting  values  in  a  reasonable  neighborhood  of  the  global  optimum. 

The  second  remark  is  that  the  ML,  and  indeed  the  MDL,  formulations  are  not  based  upon  whether  the 
underlying  polygon  is  affinely  regular,  or  even  convex.  Hence,  although  our  algorithms  may,  in  general,  yield 
non-affinely  regular  or  even  nonconvex  solutions,  they  generally  perform  best  when  the  underlying  polygon 
is  itself  affinely  regular;  this  being  a  direct  consequence  of  the  fact  that  the  initial  guess  algorithm  always 
yields  an  affinely  regular  polygon.  In  Section  3.4  we  report  studies  in  which  the  Initial  Guess  algorithm 
is  used  to  produce  a  starting  guess  for  the  optimization  routines.  In  these  studies,  we  show  that  although 
the  performance  of  the  overall  algorithm  does  degrade  somewhat,  this  degradation  in  performance  is  not 
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significant. 

In  order  to  quantify  some  measure  of  performance  of  our  proposed  reconstruction  algorithms,  we  first 
need  to  define  an  appropriate  notion  of  signal-to-noise  ratio  (SNR).  We  define  the  SNR  per  sample  as 


SNR  =  10  log10 


E  i,j92(ti,0j)/d 


(40) 


where  d  =  m  x  n  is  the  total  number  of  observations,  and  er2  is  the  variance  of  the  i.i.d.  noise  w  in  the 
projection  observations  (2). 

In  all  our  simulations  the  reconstruction  error  is  measured  in  terms  of  the  percent  Hausdorff  distance  [32] 
between  the  estimate  and  the  true  polygon  or  shape.  The  Hausdorff  metric  is  a  proper  notion  of  “distance” 
between  two  nonempty  compact  sets  and  it  is  defined  as  follows.  Let  d(p*,S)  denote  the  minimum  distance 
between  the  point  pi  and  the  compact  set  S: 


d(p*,S)  =  inf{||p*  —  p||  \p&S}. 
Define  the  e-neighborhood  of  the  set  S  as 

S(e)  =  {p  |  d(p,S)  <  e}. 


(41) 


(42) 


Now  given  two  non-empty  compact  sets,  Si  and  S2,  the  Hausdorff  distance  between  them  is  defined  as: 


H(S!,S2)  =  inf{e  |  Si  C  S(2e)  and  S2  C  s[e)} 


(43) 


In  essence,  the  Hausdorff  metric  is  a  measure  of  the  largest  distance  by  which  the  sets  Si  and  S2  differ.  The 
percent  Hausdorff  distance  between  the  true  object  S  and  the  reconstruction  S  is  now  defined  as 


Percent  Error  =  100%  x 


H{S,  S) 
H(0,S) 


(44) 


where  O  denotes  the  set  composed  of  the  single  point  at  the  origin,  so  that  if  S  contains  the  origin,  H(0 ,  S) 
is  the  maximal  distance  of  a  point  in  the  set  to  the  origin  and  thus  a  measure  of  the  set’s  size. 

In  all  experiments  that  follow,  the  field  of  view  (the  extent  of  measurement  in  the  variable  t)  is  taken 
to  be  twice  the  maximum  width  of  the  true  object.  While  it  is  true  that  significantly  enlarging  the  field  of 
view  will  certainly  degrade  the  performance  of  the  proposed  algorithms,  our  choice  to  fix  this  at  twice  the 
size  of  the  true  object  is  not  unjustified  or  arbitrary.  In  fact,  we  are  assuming  that  a  mechanish  exists,  such 
as  that  reported  in  [7],  whereby  a  rough  estimate  of  the  maximal  support  of  the  object  may  be  estimated 
from  the  raw  projection  data  in  a  statistically  optimal  fashion. 

The  numerical  algorithm  used  to  solve  the  optimization  problems  presented  in  this  paper  was  the  Nelder- 
Mead  simplex  algorithm  [34] .  The  specific  stopping  criterion  was  to  declare  a  solution  after  the  cost  failed  to 
decrease  by  more  than  10-4  in  50  iterations,  or  after  500  iterations  of  the  Nelder-Mead  algorithm;  whichever 
came  first. 

Finally,  as  for  the  assumption  that  the  underlying  object  is  simply-connected,  we  did  not  explicitly 
enforce  this  on  the  solution.  Having  obtained  a  solution  (vertices)  from  the  numerical  search  algorithm,  we 
essentially  connected  the  vertices  in  such  a  way  as  to  minimize  the  cost. 


3.1  ML  Based  Reconstruction 

Here  we  present  examples  and  performance  analyses  of  the  ML  based  reconstruction  method  (4).  As  we 
alluded  to  before,  the  studies  reported  here  correspond  to  a  local  exploration  of  the  ML  cost  function.  To 
this  end,  we  set  the  initial  guess  equal  to  the  true  object.  While  there  is  no  guarantee  that  this  choice  of 
initial  guess  will  result  in  convergence  to  the  global  optimum  (especially  given  the  presence  of  noise  in  the 
data),  the  resulting  reconstructions  suggest  that  this  initialization  typically  yields  solutions  that  are  at  least 
very  near  the  global  optimum. 
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Figure  4:  Triangle  example.  SNR=  0  dB,  50  views,  20  samples/view:  True(-),  Reconstruction^  -).  %Error 
=  7.2 

3.1.1  Sample  Reconstructions 

In  Figures  4  and  5  we  show  optimal  reconstructions  of  a  triangle  and  a  hexagon,  respectively,  based  on 
the  ML  criterion.  The  true  polygon,  in  each  case  is  depicted  in  solid  lines  while  the  estimate  is  shown 
in  dashed  lines.  For  both  objects,  1000  noisy  projection  samples  were  collected  in  the  form  of  50  equally 
spaced  projections  in  the  interval  (0, 7r]  (m=50),  and  20  samples  per  projection  (n=20).  The  field  of  view 
(extent  of  measurements  in  the  variable  t)  was  chosen  as  twice  the  maximum  width  of  the  true  object  in 
each  case.  For  each  of  these  data  sets  the  variance  of  the  noise  in  (2)  was  set  so  that  the  SNR  given  by  (40) 
was  equal  to  0.  The  typical  behavior  of  the  optimal  ML  based  reconstructions  in  the  projection  space  can 
be  seen  in  Figure  6,  which  corresponds  to  the  hexagon  of  Figure  5.  The  top  image  of  this  figure  shows  the 
underlying  projection  function  g(ti,9j,V*)  of  (2)  for  the  hexagon,  while  the  middle  image  shows  the  noisy 
observed  data  Yhj.  The  object  is  difficult  to  distinguish  due  to  the  noise  in  the  image.  The  bottom  image 
shows  the  reconstructed  projections  corresponding  to  the  optimal  estimate  g(ti,9j,V),  which  are  virtually 
indistinguishable  from  those  corresponding  to  the  true  object.  Figure  7  shows  the  best  FBP  reconstruction 
of  the  hexagon  used  in  Figure  5  based  on  4096  projection  samples  of  the  same  SNR  (0)  (64  angles  with  64 
samples  per  angle).  For  comparison,  the  reconstruction  from  this  data  using  our  algorithm  is  shown  in  Figure 
8.  (Note  here  that  what  constitutes  the  “best”  FBP  is  somewhat  subjective.  We  tried  many  different  filters 
and  visually,  the  best  reconstruction  was  obtained  with  a  Butterworth  filter  of  order  3  with  0.15  normalized 
cutoff  frequency.) 

Note  that  the  number  of  samples  per  projection  used  in  this  reconstruction  is  actually  more  than  the 
number  used  to  produce  the  ML-based  reconstruction  in  Figure  5.  The  increase  in  sampling  was  necessary 
because  CBP  produces  severe  artifacts  if  the  number  of  views  exceeds  the  number  of  samples  per  view  [5] . 
The  ML  approach  has  no  such  difficulties,  as  we  will  see  in  the  next  section  where  we  examine  performance. 

In  contrast  to  the  ML-based  reconstruction,  the  details  of  the  hexagon  are  corrupted  in  the  FBP  recon¬ 
struction.  In  addition,  there  are  spurious  features  in  the  FBP  reconstructions  and  perhaps  most  importantly, 
to  extract  a  binary  object  from  the  FBP  reconstruction,  we  would  need  to  threshold  the  image  or  perform 
edge  detection  on  it.  Neither  of  these  postprocessing  steps  are  easily  interpretable  in  an  optimal  estimation 
framework  and,  of  course,  they  incur  even  more  computational  costs. 
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Figure  5:  Hexagon  example.  SNR=  0  dB  50  views,  20  samples/view:  True(-),  Reconstruction^  -). 
%Error=9.6 


3.1.2  Effect  of  Noise  on  Performance 

The  average  performance  of  the  ML  based  reconstructions  is  presented  through  several  Monte-Carlo  studies. 
Again,  for  these  experiments  an  initial  condition  equal  to  the  true  object  was  used  in  each  case  to  ensure  us 
of  obtaining  the  actual  ML  estimates.  The  first  study  establishes  average  reconstruction  errors  at  various 
SNR’s  for  a  fixed  number  of  data  samples.  The  purpose  of  these  simulations  is  to  demonstrate  that  the 
ML-based  reconstructions  are  robust  to  the  quality  of  the  data  used  for  a  wide  range  of  SNR’s.  The  same 
two  polygons  as  in  Figures  4  and  5  were  chosen  as  the  underlying  objects.  Again,  in  each  case,  1000  samples 
of  the  projections  of  these  objects  were  collected  in  the  form  of  50  equally  spaced  projections  in  the  interval 
(0, 7r]  (m=50),  and  20  samples  per  projection  (n=20)  while  the  field  of  view  (extent  of  measurements  in  the 
variable  t)  was  chosen  as  twice  the  maximum  width  of  the  object  in  each  case.  The  samples  g(ti,  9j ,  V*)  were 
then  corrupted  by  Gaussian  white  noise  w(ti,6j)  of  different  intensities  to  yield  data  sets  at  several  SNR’s. 
At  each  SNR,  100  reconstructions  were  done  using  independent  sample  paths  of  the  corrupting  white  noise. 
The  average  reconstruction  error  was  then  computed  and  is  displayed  against  the  SNR  in  Figure  9.  The 
error  bars  denote  the  95%  confidence  intervals  for  the  computed  mean  values. 

The  percent  error  in  these  reconstructions  increases  with  decreasing  SNR,  as  one  would  expect.  In  fact, 
the  graph  shows  that,  at  least  in  the  examined  SNR  range  of  —4.35  to  +4.35  dB,  the  relation  between 
percent  error  and  SNR  is  roughly  linear  in  the  cases  of  the  triangle  and  the  hexagon.  This  suggests  that 
given  a  good  initial  guess,  the  performance  of  our  algorithm  does  not  degrade  very  fast  with  decaying  SNR, 
demonstrating  the  robustness  to  noise  of  such  object-based  optimal  ML  estimates. 

3.1.3  Effect  of  Sampling  on  Performance 

Here  the  performance  of  our  ML-based  estimates  with  respect  to  both  the  number  of  available  data  samples 
and  their  distribution  is  studied.  One  would  naturally  expect  that  as  the  number  of  available  data  points 
decreases,  the  reconstruction  error  should  increase.  The  main  aim  of  these  simulations  is  to  demonstrate 
that  given  a  good  initial  guess,  the  ML-based  reconstructions  are  robust  to  both  the  quantity  and  the 
distribution  of  data  over  a  wide  range  of  SNR’s.  In  particular,  reasonable  estimates  are  produced  even  with 
a  drastic  reduction  of  data  and,  unlike  the  behavior  seen  in  FBP  reconstructions,  the  ML  estimates  display 
no  catastrophic  degradation  as  the  samples  per  angular  view  are  reduced.  The  relative  sensitivity  of  the  ML 
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Figure  6:  From  Top  to  Bottom:  Sinograms  with  50  projections  and  20  samples  per  projection  of  I)  Noiseless 
Hexagon,  II)  Noisy  data  at  0  dB  III)  Reconstructed  Hexagon.  In  each  of  these  images,  the  horizontal  axis  is 
9,  the  vertical  axis  is  t,  and  the  intensity  values  are  the  values  of  the  corresponding  projections  mapped  to 
the  grayscale  range  of  [0,255] 


estimates  to  density  of  samples  of  g(t ,  9 ,  V*)  in  t  and  9  is  also  discussed,  providing  information  of  use  for  the 
design  of  sampling  strategies. 

The  true  hexagon  used  in  Figure  5  was  again  used  as  the  underlying  object.  As  before,  an  initial  condition 
equal  to  the  true  object  was  used  for  each  of  experiments.  A  series  of  Monte-Carlo  simulations  (50  runs  for 
each  sampling  configuration)  were  then  performed  at  various  SNR’s  to  observe  the  effect  of  sparse  projections 
and  sparse  sampling  in  each  projection.  In  Figure  10,  the  percent  Hausdorff  reconstruction  error  is  plotted 
versus  the  number  of  angular  views  for  SNR’s  of  0,  4.35,  and  8.7  dB,  while  the  number  of  samples  per  view 
was  fixed  at  50.  With  a  modest  50  samples  per  view,  all  three  curves  fall  below  10%  reconstruction  error 
when  the  number  of  views  is  greater  than  about  10.  This  is  only  500  total  observations,  many  of  which  do 
not  contain  the  object  at  all  (since  the  field  of  view  is  twice  as  large  as  the  object).  Furthermore,  as  the 
number  of  angular  views  is  decreased  from  100  to  10,  only  a  marginal  increase  in  the  reconstruction  error  is 
observed.  These  observations  testify  to  the  robustness  of  the  ML  algorithm  with  respect  to  the  number  of 
views  when  a  good  initial  guess  is  given. 

In  Figure  11,  the  dual  case  is  presented.  In  this  figure  the  percent  Hausdorff  reconstruction  error  is 
plotted  versus  the  number  of  samples  per  view  for  SNR’s  of  0,  4.35,  and  8.7  dB,  while  the  number  of  angular 
views  was  fixed  at  50.  With  50  angular  views,  all  curves  fall  below  10%  reconstruction  error  when  the 
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Figure  7:  Sample  reconstruction  of  a  hexagon  at  0  dB  SNR  using  FBP:  64  views,  64  samples  per  view 


Figure  8:  Sample  reconstruction  of  a  Hexagon  at  0  dB  SNR  64  views,  64  samples/view:  True(-), 
Reconstruction (-  -) 


number  of  samples  per  view  is  greater  than  only  10.  Also,  as  the  number  of  samples  per  view  is  decreased 
from  100  to  10,  again  only  a  marginal  increase  in  the  reconstruction  error  is  observed.  This  behavior  shows 
that  the  ML  algorithm  is  robust  with  respect  to  the  number  of  samples  per  view  when  a  good  initial  guess 
is  given.  Note  that  for  a  fixed  sampling  strategy,  the  reconstruction  error  increases  only  slightly  as  the  SNR 
is  decreased  over  a  wide  range.  For  instance,  in  Figure  10,  with  40  angular  views  and  50  samples  per  view, 
the  percent  error  is  reduced  only  about  5%  while  the  SNR  goes  from  0  to  8.7  dB. 

Finally,  it  is  noteworthy  that  the  reconstruction  error  enjoys  a  dramatic  improvement  for  all  SNR’s  (0, 
4.35,  and  8.7  dB)  when  the  number  of  samples  per  view  is  increased  from  5  to  10.  This  improvement  is 
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Figure  9:  Mean  performance  curves  for  ML  reconstructions  of  a  triangle  and  a  hexagon 


Figure  10:  Performance  as  a  function  of  Number  of  Views 


more  significant  than  that  observed  in  Figure  10  when  the  number  of  views  is  increased  from  5  to  10.  This 
behavior  indicates  that  in  a  scenario  where  only  a  small  (fixed)  number  of  sample  points  can  be  collected,  it 
is  more  beneficial  to  have  more  samples  per  view  rather  than  more  views. 

3.2  MDL  Reconstructions 

Here  we  will  examine  reconstruction  under  the  MDL  criterion  of  (7)  where  we  now  assume  that  the  number 
of  sides  of  the  reconstructed  polygon  is  unknown.  In  particular,  the  reconstruction  experiments  for  the 
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Figure  11:  Performance  as  a  Function  of  Number  of  Samples  per  View 


hexagon  in  Figure  5  were  repeated  at  SNR=0  dB  assuming  no  knowledge  of  the  number  of  sides.  The  MDL 
criterion  was  employed  to  estimate  the  optimal  number  of  sides.  As  in  the  ML  algorithm,  it  is  important  to 
find  a  good  initial  guess  for  the  MDL  algorithm  as  well.  The  problem  is  twofold.  First,  a  reasonable  guess 
must  be  made  as  to  the  appropriate  range  of  the  number  of  sides.  We  picked  a  fairly  small  range  for  the 
number  of  sides  of  the  reconstruction;  typically,  3  to  10  sides.  Next,  for  each  number  of  sides,  the  Initial 
Guess  algorithm  was  used  to  produce  an  initial  guess  to  the  optimization  routine.  The  method  for  selecting 
the  range  of  the  number  of  sides  is  ad  hoc,  but  was  shown  to  be  reliable  in  the  sense  that  for  our  simulations, 
the  MDL  cost  never  showed  local  or  global  minima  for  convex  objects  with  number  of  sides  larger  than  10. 
Figure  12  shows  a  plot  of  the  MDL  cost  corresponding  to  the  expression  in  (7)  versus  the  number  of  sides 
for  a  sample  reconstruction  of  the  hexagon  in  Figure  5.  It  can  be  seen  that  the  minimum  occurs  at  N  =  6, 
demonstrating  that  the  optimal  MDL  reconstruction  will  consist  of  6-sides.  Indeed  this  number  coincides 
with  the  true  number  of  sides  of  the  underlying  object.  The  optimal  MDL  estimate  is  thus  exactly  the 
optimal  ML  estimate  for  this  data  set  presented  before. 

3.3  Polygonal  Reconstruction  of  Non-polygonal  Objects 

In  this  section  we  wish  to  study  the  robustness  of  MDL-based  estimates  when  the  underlying,  true  object  is 
non-polygonal.  First  we  examine  the  case  of  an  elliptical  object.  We  use  the  MDL  formulation  presented  in 
the  previous  section  and  study  the  behavior  of  the  optimal  reconstructions  at  two  different  SNR’s.  To  this 
end,  let  the  true  object  (that  which  generated  the  data)  be  a  binary  ellipse  whose  boundary  is  given  by: 

+  =  (45) 

The  above  relation  defines  an  ellipse  centered  at  the  point  (1/2,— 1/2)  whose  major  and  minor  axes  are 
aligned  with  the  coordinate  axes  with  lengths  1  and  3/2,  respectively. 

One  thousand  (1000)  noisy  samples  of  the  Radon  transform  of  this  ellipse  were  generated  (m=50  equally 
spaced  angular  views  in  (0,7r],  and  n=20  samples  per  view)  at  SNR’s  of  0  and  2.17  dB  respectively  for  50 
different  sample  paths  of  the  corrupting  noise.  For  each  set  of  data,  reconstructions  were  performed  using 
the  ML  algorithm  with  3,  4,  5,  6,  7,  and  8  sides  together  with  the  initial  guess  algorithm.  The  MDL  cost  in 
(7)  was  then  computed  for  each  of  these  reconstructions.  The  ensemble  mean  of  this  cost  over  the  50  runs, 
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Figure  12:  Cost  vs  number  of  sides  for  the  hexagon  in  Figure  5 


5-sided  reconstruction,  SNR=2.17  dB 
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6-sided  reconstruction,  SNR=0  dB 


Figure  13:  Minimum  MDL  costs  and  Sample  Reconstructions  for  an  Ellipse 


for  each  value  of  N,  is  presented  in  the  top  part  of  Figure  13.  The  error  bars  denote  the  95%  confidence 
intervals  for  the  computed  mean  values.  The  top  left  curve  corresponding  to  the  SNR=  2.17  dB  case  displays 
its  minimum  at  N  =  5.  This  behavior  indicates  that  the  average  optimal  MDL  reconstruction  uses  5  sides  at 
this  noise  level.  A  corresponding  typical  such  five-sided  reconstruction  of  the  ellipse  is  displayed  on  the  lower 
left  plot  of  Figure  13  together  with  the  true  ellipse.  The  upper  right  curve  corresponding  to  the  SNR=  0  dB 
case  displays  its  minimum  at  N  =  6  which  indicates  that  the  average  optimal  MDL  reconstruction  for  this 
case  uses  6  sides.  The  MDL  cost  curve  for  this  lower  SNR  case  has  now  become  quite  flat  however,  showing 
that  the  reconstruction  with  N  from  4  to  6  are  all  about  equally  explanatory  of  the  data.  Although  the  curves 
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9-sided 


Figure  14:  True  object  (-),  reconstruction  initial  guess  (o)  picked  using  Initial  Guess  Algorithm,  SNR= 
4.35  dB. 


for  both  cases  demonstrate  the  ability  of  the  MDL  procedure  to  capture  the  shape’s  complexity  through  its 
choice  of  N,  this  behavior  suggests  that  with  increasing  noise  intensity,  an  MDL-based  estimate  becomes  less 
sensitive  to  the  precise  level  of  complexity  of  the  reconstruction,  as  we  would  expect.  Apparently,  in  high 
noise  situations  the  differences  between  these  reconstructions  that  would  be  apparent  in  high  SNR  scenarios 
are  masked.  As  the  noise  level  increases,  these  fine  distinctions  are  unimportant  or  not  supported  by  the 
data.  A  typical  6-sided  reconstruction  is  also  displayed  in  the  lower  right  plot  of  Figure  13  along  with  the 
true  ellipse. 

As  another  example  of  the  use  of  the  proposed  algorithms,  we  choose  an  object  that  is  non-polygonal  and 
non-convex.  In  Figure  14,  reconstructions  of  this  object  are  shown  based  on  20  equally  spaced  projections 
with  50  samples  per  projection,  at  a  signal  to  noise  ratio  of  4.35  dB.  In  these  reconstructions,  a  slight 
variant  of  the  initial  guess  algorithm  was  used  to  generate  the  starting  polygons  as  we  will  describe  shortly. 
Figure  15,  meanwhile,  contains  the  reconstruction  produced  by  FBP  using  the  same  data  set.  As  can  be 
seen,  when  using  seven  or  more  sides,  the  underlying  object  can  be  captured  more  accurately  and  without 
spurious  features  through  the  use  of  our  algorithm. 

To  pick  the  “best”  number  of  sides,  in  Figure  16  we  show  the  MDL  cost  curve  for  the  reconstructions  of 
the  kidney-shaped  object  at  SNR  =  4.35  dB  shown  in  Figure  14.  In  applying  our  algorithm  to  this  non-convex 
object,  we  found  that  we  can  improve  the  resulting  reconstructions  by  slightly  perturbing  the  initialization 
produced  by  the  proposed  Initial  Guess  Algorithm  by  rotations  through  small  angles.  In  particular,  recall 
that  according  to  Result  1,  any  orthogonal  transformation  of  the  initial  polygon  Unit  produced  by  the  Initial 
Guess  Algorithm  is  a  “valid”  starting  guess  in  the  sense  that  it  will  have  the  same  first  three  moments  as 
Unit-  Taking  advantage  of  this  property,  for  a  given  number  of  sides,  we  applied  several  small  rotations  to 
the  Unit  produced  by  the  Initial  Guess  Algorithm,  and  carried  out  the  ML  optimization  problem  with  these 
resulting  initializations.  We  then  picked  the  solutions  which  resulted  in  the  lowest  ML  cost.  Using  this 
procedure  for  each  number  of  sides,  we  calculated  the  MDL  cost  values  shown  in  Figure  16.  These  “best” 
reconstructions  are  shown  in  Figure  14. 2  We  note  that,  according  to  Figure  16,  the  minimum  of  the  MDL 
cost  corresponds  to  a  reconstruction  with  7  sides.  As  is  apparent  in  Figure  14,  while  using  a  higher  number 

2The  process  of  applying  an  orthogonal  transformation  to  the  output  of  the  Initial  Guess  Algorithm,  before  using  this  in 
the  reconstruction  algorithm,  is  one  that  can  perhaps  be  applied  in  general.  However,  we  found  this  modification  to  be  needed 
only  when  rather  complex,  non-convex  and  non-polygonal  objects  were  being  reconstructed. 
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Figure  15:  FBP  Reconstruction  of  non-polygonal,  non-convex  Object:  3rd  order  Butterworth  filter  with  0.15 
normalized  cutoff  frequency,  SNR=4.35  dB 


Figure  16:  MDL  cost  curve  for  the  reconstruction  of  the  kidney-shaped  object  with  the  Initial  Guess  Algo¬ 
rithm  used. 


of  vertices  beyond  7  does  improve  the  reconstructions  somewhat,  it  does  not  yield  significantly  better  results. 
This  fact  is  directly  reflected  in  the  MDL  cost  curve  becoming  increasingly  flat  beyond  7  sides. 


24 


350 


300 


.  250 


o  200 
-a 


a 

X 


;  150 


■  100 


50 


Outlier  — > 


Outlier  — > 


20  40  60  80  100 

Realizations  Number 


Figure  17:  A  sample  path  of  the  reconstruction  error  at  SNR=0  dB 

3.4  Initial  Guess  Algorithm 

In  this  section  we  present  some  sample  reconstructions  and  performance  plots  for  which  we  use  the  Initial 
Guess  algorithm  for  generating  a  starting  point  to  the  nonlinear  optimization  (3).  These  experiments  aim 
to  show  that  the  Initial  Guess  Algorithm  does  indeed  provide  us  with  a  starting  guess  that  in  the  great 
majority  of  the  cases  is  near  the  global  optimum.  This  will,  of  course,  not  guarantee  convergence  to  the 
global  optimum  solution,  but  as  we  shall  see,  the  use  of  the  Initial  Guess  Algorithm  does  at  least  result  in 
convergence  to  a  local  minimum  of  the  cost  which  is  almost  always  quite  near  the  global  optimum. 

To  study  the  average  performance  of  the  ML  algorithm  using  the  output  of  the  Initial  Guess  Algorithm, 
a  Monte-Carlo  simulation  was  done  for  the  reconstruction  of  the  hexagon  shown  in  Figure  5.  100  recon¬ 
structions  were  carried  out  for  different  realizations  of  the  noise  at  various  SNR’s,  each  with  1000  projection 
samples  as  before  (50  projections  and  20  samples  per  projection).  For  each  SNR,  on  average  less  than  5 
percent  of  the  reconstructions  (i.e.  5  out  of  100  sample  reconstructions)  had  very  large  reconstruction  errors 
(we  call  these  instances  outliers).  Figure  17  shows  the  reconstruction  errors  for  100  realizations  of  the  noise 
at  0  dB.  The  outliers  are  clearly  visible. 

Figure  17  indicates  that  in  a  few  instances,  the  reconstructions  were  essentially  at  local  minima  very  far 
from  the  global  minimum  of  the  cost.  In  our  experience,  these  outliers  occur  most  frequently  when  poor 
estimates  of  the  moments  of  order  2  are  obtained  from  the  noisy  data.  Note  that  the  second  order  moments 
are  used  in  the  Initial  Guess  Algorithm  only  if  the  corresponding  inertia  matrix  obtained  from  them  is  strictly 
positive  definite.  The  Initial  Guess  Algorithm  decides  whether  to  use  the  second  order  moments  or  not  solely 
on  the  basis  of  this  positive  definiteness  and  regardless  of  how  close  the  inertia  matrix  may  be  to  negative-  or 
indefiniteness.  Hence,  the  outliers  occur  in  those  rare  instances  when  the  estimated  inertia  matrix  happens 
to  be  a  very  poor  estimate,  but  yet  positive  definite  (and  hence  used  in  the  Initial  Guess  Algorithm).  This 
phenomenon,  in  turn,  seems  to  occur  when  relatively  few  samples  per  projection  are  available. 

Figure  18  shows  the  mean  percent  error  in  the  Monte-Carlo  runs  after  the  removal  of  the  outliers.  The 
outliers  were  removed  from  the  ensemble  and  the  results  of  the  remaining  realizations  were  averaged  to  yield 
the  values  in  Figure  18.  That  is  to  say  that  if  3  out  of  100  realizations  led  to  outliers,  then  only  those  97 
results  which  seemed  “reasonable”  were  used  in  computing  the  ensemble  average.  Whether  the  result  of  a 
run  was  deemed  reasonable  or  not  was  decided  by  comparing  the  resulting  percent  error  to  the  ensemble 
median  reconstruction  error  for  all  100  runs.  In  particular,  if  the  percent  error  for  a  run  was  larger  than 
one  standard  deviation  away  from  the  median,  the  run  was  declared  an  outlier.  In  the  case  of  Figure  17, 
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Figure  18:  Percent  error  for  Hexagon  vs  SNR  after  outlier  removal 

using  the  computed  median  value  of  17.2,  and  standard  deviation  of  33.2  a  threshold  level  of  50.4  was  chosen 
above  which  outliers  were  declared. 

The  resulting  “mean”  performance  is  plotted  here  to  show  the  average  performance  without  the  effect  of 
the  outliers.  It  can  be  seen,  upon  comparing  Figure  18  with  the  corresponding  performance  curve  in  Figure 
9  that  the  performance  of  the  ML  algorithm  using  the  output  of  the  Initial  Guess  algorithm  still  suffers 
even  after  discounting  the  obvious  outliers.  This  means  that  instances  of  convergence  to  local  minima  still 
occur,  but  note  that,  at  least  from  a  visual  standpoint,  the  average  performance  after  the  removal  of  outliers 
is  not  significantly  different  from  the  average  performance  with  the  actual  polygon  as  the  initial  guess.  In 
particular,  the  degradation  in  performance  here  is  roughly  7  percentage  points  in  the  Hausdorff  norm  over  the 
given  SNR  range.  This  corresponds  to  a  small  visual  error  as  can  be  seen  in  Figure  4.  From  this  observation, 
we  conclude  that  even  though  the  initial  guess  algorithm  does  not  always  lead  to  convergence  to  the  global 
minimum,  it  almost  always  leads  to,  at  least,  a  local  minimum  that  is  fairly  close  to  the  global  minimum  of 
the  cost  function.  Typical  reconstruction  at  local  minima  which  are  close  to  the  global  minimum  of  the  cost 
are  shown  in  Figure  19  for  SNR=  0  dB  and  in  Figure  20  for  SNR=  4.35  dB. 

4  Conclusions 

In  this  paper,  we  studied  statistical  techniques  for  the  reconstruction  of  binary  polygonal  objects.  In  partic¬ 
ular,  we  focused  on  the  reconstruction  of  binary  polygonal  objects.  The  reconstruction  of  such  objects  was 
posed  as  a  parameter  estimation  problem  for  which  the  Maximum  Likelihood  technique  was  proposed  as  a 
solution.  In  contrast  to  the  classical  techniques,  such  as  FBP,  the  ML  based  reconstructions  showed  great 
robustness  to  noise  and  data  loss  and  distribution  when  a  good  initial  guess  was  available.  The  drawback 
of  such  ML-based  formulations  is  that  the  resulting  optimization  problems  are  highly  non-linear  and  thus  a 
good  initial  guess  is  necessary  to  ensure  convergence  of  optimization  routines  to  a  solution  that  is  at  least 
near  the  true  ML  estimate.  To  this  end,  an  algorithm  was  presented  for  computing  such  a  reasonable  initial 
guess  using  moments  of  the  object  which  are  estimated  directly  from  the  projection  data.  While  estimation 
of  a  function  from  its  moments  is,  in  general,  a  difficult  and  ill-posed  problem,  we  avoid  these  problems  by 
using  the  noisy  estimated  moments  only  to  guide  a  coarse  object  estimate.  This  estimate,  in  turn,  mitigates 
the  difficulties  associated  with  the  non-linearities  of  the  optimal  ML  statistical  approach.  The  efficacy  of 
this  moment  based  initial  guess  algorithm  was  demonstrated  over  a  range  of  SNR’s. 
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True(-),  Reconstruction^  .),  Initial  Guess(+) 


Figure  19:  A  typical  reconstruction  at  a  local  minimum  with  SNR=  0  dB 


True(-),  Reconstruction^  .),  Initial  Guess(+) 


Figure  20:  A  typical  reconstruction  at  a  local  minimum  with  SNR=  4.35  dB 


If  the  number  of  parameters  describing  the  underlying  object  are  not  known,  a  Minimum  Description 
Length  criterion  can  be  employed  that  simply  generalizes  the  ML  framework  to  penalize  the  use  of  an 
excessively  large  number  of  parameters  for  the  reconstruction.  The  MDL  approach  was  shown  to  work 
successfully  in  estimating  the  number  of  sides  and  the  underlying  object  itself  for  low  signal-to-noise  ratio 
situations  and  for  a  variety  of  sampling  scenarios.  It  was  further  demonstrated  that  if  the  underlying  object 
is  not  polygonal,  but  still  binary,  the  proposed  ML  and  MDL  algorithms  are  still  capable  of  producing 
polygonal  reconstructions  which  reasonably  capture  the  object  shape  in  the  presence  of  high  noise  intensity. 

In  this  work  we  have  focused  on  the  reconstruction  of  binary  polygonal  objects  parameterized  by  their 
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vertices.  The  ML  and  MDL-based  techniques  used  here  may  also  be  applied  to  more  general  object  pa- 
rameterizations.  In  particular,  while  we  used  the  (estimated)  moments  of  the  object  only  as  the  basis  for 
generating  an  initial  guess,  it  is,  in  some  cases,  possible  to  actually  parameterize  the  object  entirely  through 
its  moments.  For  instance,  Davis  [35]  has  shown  that  a  triangle  in  the  plane  is  uniquely  determined  by  its 
moments  up  to  order  3,  while  in  [27,  33]  we  have  generalized  this  result  to  show  that  the  vertices  of  any 
simply  connected  nondegenerate  TV-gon  are  uniquely  determined  by  its  moments  up  to  order  2 N  —  3. 

More  generally,  a  square  integrable  function  defined  over  a  compact  region  of  the  plane  is  completely 
determined  by  the  entire  set  of  its  moments  [31,  29,  28].  In  reality  we  will  only  have  access  to  a  finite 
set  of  these  moments  and  these  numbers,  coming  from  estimates,  will  themselves  be  inexact  and  noisy. 
While  estimation  of  the  moments  of  a  function  based  on  its  projections  is  a  convenient  linear  problem, 
inversion  of  the  resulting  finite  set  of  moments  to  obtain  the  underlying  function  estimate  is  a  difficult 
and  ill-posed  problem.  These  observations  suggest  a  spectrum  of  ways  in  which  to  use  moments  in  our 
reconstruction  problems.  At  one  extreme,  only  a  few  moments  are  used  in  a  sub-optimal  way  to  generate 
a  simple  initialization  for  solution  of  a  hard,  non-linear  estimation  problem.  At  the  other  extreme,  the 
moments  are  themselves  used  in  an  optimal  reconstruction  scheme.  In  [27,  26]  we  have  studied  regularized 
variational  formulations  for  the  reconstruction  of  a  square  integrable  function  from  noisy  estimates  of  a  finite 
number  of  its  moments.  We  have  also  studied  [33]  array-processing-based  algorithms  for  the  reconstruction 
of  binary  polygonal  objects  from  a  finite  number  of  their  moments. 
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A  Theoretical  Results  on  the  Initial  Guess  Algorithm 

In  this  section  we  present  some  theoretical  justification  for  the  initial  guess  algorithm.  To  start,  we  state 
some  elementary  properties  of  unit  area  polygons  Vle{  (N)  whose  vertices  are  the  scaled  Nth  roots  of  unity 
(in  counter-clockwise  direction)  as  defined  by  (18).  From  [27],  it  is  a  matter  of  some  algebraic  manipulations 
to  show  that  the  regular  polygon  Vref  (N)  has  moments  of  up  to  order  2  given  by 


Moo(Kef  (AO)  =  1  (46) 

MloTref(A0)=  MOlTref(AO)  =0  (47) 

M2o(Kef  (N))=  M02(Kef  (N))  =  ^  ^  =  kN  (48) 

Mil  (Kef  (N))  =  0  (49) 

Now  let  VJnit  be  an  affine  transformation  of  Vle{  as 

Knit  =  TKef(iV)  +  [C\C\  ■■■  \C\  (50) 


for  some  linear  transformation  L  and  some  2x1  vector  C.  Let  O(Knit)  denote  the  closed,  binary  polygonal 
region  enclosed  by  the  7V-gon  Vjn;t.  Now  by  considering  the  change  of  variables  z  =  Lu,  and  dropping  the 
explicit  dependences  on  N  we  have 


Moo  (Knit) 


/  [  dz, 

■1  JO(V: init) 

(51) 

f  f  det(L)  du, 

./  JO(Vte{) 

(52) 

Moo  (Kef)  det(L)  =  |  det(L) 

(53) 
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Similarly,  we  get 


[Mio(Mnit)  M0i(Vinit)]T  =  (L[/i lo(Vref)  Moi(Kef)]T  +  C)  |  det(X)|  =  I  det(L)|C  (54) 


and 

J(Mnit)  =  ( Ll(Viel)LT  +  CCT)\  det(i)|  =  (/cjvLLt  +  CCT)|  det(L)|  (55) 


where  for  any  N-gon  V  we  write 


1{V)  = 


V20(V)  /in(F) 

m(V)  ^{V) 


(56) 


This  proves  relations  (21),  (22),  and  (23).  We  next  establish  an  explicit  description  of  the  set  of  all  affinely 
regular  IV-gons  with  a  fixed  set  of  moments  up  to  order  2.  In  order  to  do  this,  we  first  need  to  prove  a 
lemma. 


Lemma  1  For  every  N-gon  V  with  moments  /zoo,  Mio  =  0,  /zoi  =  0,  /z 20,  /in,  M02,  such  that  the  inertia 
matrix  X  satisfies  det(X)  =  fc^/z g0,  there  exists  a  matrix  L,  unique  up  to  some  orthogonal  transformation, 
such  that  V  =  LV, ref. 


Proof:  The  assumptions  that  /iio  =  0  and  /ioi  =  0  are  made  without  loss  of  generality  and  to  facilitate 

the  presentation  of  the  proof.  Having  said  this,  we  define  L  as  the  scaled  (unique)  square  root  of  X  as  follows. 
First,  write  the  following  eigendecomposition 


X 


US2UT, 


where  JJ  is  orthogonal  and  S  has  unit  determinant.  Define  L  as 


(57) 


L  =  VJ^US.  (58) 

The  moments  of  V  =  LVref  are  then  given  by 

Poo(V)  =  Moo  (59) 

Mio(V)  =  Moiiy)  =  0  (60) 

X(V)  =  fc/vXXT|  det(X)|  (61) 

Note  that 

det(X)  =  fcjv/r  00  (62) 


as  required.  If  L  is  replaced  by  LT  where  T  is  any  2x2  orthogonal  transformation,  the  same  moments  are 
obtained.  Hence  the  lemma  is  established.  □ 


Given  this  lemma,  we  obtain  an  interesting  geometric  representation  of  all  affinely  regular  IV-gons  that 
have  a  prespecified  set  of  moments  of  up  to  order  2.  This  characterization  is  given  by  Result  1,  on  page  13, 
which  we  prove  next. 


Proof  of  Result  1:  For  the  sake  of  simplicity,  and  without  loss  of  generality,  we  carry  out  the  proof  for 
the  case  where  all  polygons  are  centered  at  the  origin. 

Let  Si  denote  the  set  of  all  N- gons  whose  first  three  moment  sets  are  /too,  M10  =  Moi  =  0,  M20,  Mil,  Mo2- 
Let  S2  denote  the  set  of  all  IV-gons  with  vertices  on  the  ellipse  ,ltEq1'z  =  1,  and  sides  tangent  (at  their 
mid-points)  to  the  ellipse  zTEf1z  =  1.  We  show  that  Si  =  S2. 

First  consider  an  IV-gon  V  €  Si.  V  has  moments  /too,  0,  0,  X  and  therefore,  by  Lemma  1,  there  exists  an 
L  given  by  (57)  and  (58),  unique  up  to  some  orthogonal  matrix  T\  such  that  we  can  write 

V  =  (LTi)Vref(N)  (63) 
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Let  us  denote  the  IV-gons  V  and  Vref(N)  explicitly  in  terms  of  their  columns  as 


V  =  [vi  |  v2  |  •  •  •  |  vN] 
Vref(N)  =  [W!  |  W2  |  '  '  •  |  WN } 

so  that 

Vj  =  LTiWj. 

It  is  easy  to  show  from  the  definition  of  Vief(N)  that 

wj  Wj  =  Ot  n  = 


Y  sin(27r/TV)  ’ 


(Wj+ 1  +  Wj)  (wj+ 1  +  Wj)  =  A(3n  - 


Artan(7r/TV) 


Now  to  show  that  V  €  S2,  we  prove  that 

(Vj+ 1  +  Vj)T  ri„1  (Vj+ 1  +  Vj) 

2  2 


vjEo\  =  1 


=  1 


for  j  =  1,  2  •  •  •  N,  where  by  convention,  N  +  1  =  1.  Using  (57)  and  (58),  we  can  write 


vjE0\  = 


l-lOO^N  T- 7 — 1 


Vjl  Vj. 


OtN 
/ioofcv  1 

OtN  fc/VMoo 


Moo  wj  Tl  SUtUS-2UtUST1wj  . 


1 

OtN 

1. 


Similarly, 


(i’j+ 1  +  ^')T  F-i  fa'+i  +  *b)  _  lMoofcw _ 1 


-2 - ^ 


Moofe+i  +vj)  2  (vj+ 1  +  Uj) 


4/3jv 

=  1 

Hence,  V  <E  52. 

Now  assume  that  V  €  S2.  Then,  by  assumption, 


4  /3jv  fcjv  M  oo 

K'+i +  wj)TK+i +  wj) 


vjE01vj  =  1 

(vj+i  +  wj)T  *,-1  (vf+i  +  vj) 

2  Ef  2  —  1 


Define  the  vertices  of  the  related  IV- gon  Z  as 


S~  U  Vj. 

y/OlN^  00 


(64) 

(65) 

(66) 

(67) 

(68) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 


(75) 

(76) 

(77) 

(78) 

(79) 


where  S  and  U  are  given  by  the  normalized  eigendecomposition  of  I  given  in  (57).  Writing  (77)  and  (78)  in 
terms  of  Zj,  after  some  algebraic  manipulations,  we  get 


Zj  Zj  —  1 


{Zj+I  +  Zj)T  (Zj+1  +  Zj) 
2  2 
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7 r  . 


=  C0S^N] 


(80) 

(81) 


From  these  identities,  again  with  some  algebraic  manipulation,  it  easily  follows  that  Z  is  equilateral.  Specif¬ 
ically, 


lkj+i  -*jl!1/2  =  2sin(^0- 


(82) 


Since  the  above  identities  show  that  the  IV-gon  Z  is  a  regular  IV-gon  inscribed  in  the  unit  circle,  then  it 
must  be  related  to  Kef  (N)  through  a  scaling  and  some  orthogonal  transformation  T2.  In  particular, 


Z  = 


N  2-7T 

Ysin(-)T2Kef  (N) 


This  in  turn  shows  that 


or  after  solving  for  V  and  simplifying, 


Letting  L  =  ^Moo  US,  we  obtain 


V  =  yjI^UST2Vie{(N). 


V  =  (LT2)Vref(N) . 

This  last  identity  implies  that  V  has  moments  moo(K)  =  Moo ,  M 10 (K)  =  p 01  (K)  =  0,  and 

X(K)  =  kNLLT\ det(X)| 

with  det(X(K))  =  Hence  V  £  Si  and  the  result  is  established. 


(83) 

(84) 

(85) 

(86) 

(87) 

□ 


If  X  is  not  the  inertia  matrix  of  an  affinely  regular  7V-gon  then  the  L  constructed  in  the  Initial  Guess 
Algorithm  will  not  have  have  the  prescribed  inertia  matrix.  We  are,  however,  able  to  explicitly  compute  the 
approximation  error  in  the  following  way. 


Result  2  Suppose  that  the  moments  /ioo,  n  10  =  M01  =  0, 

1  = 

are  given  such  that  det(X)  =  k^p, g0  +  e  >  0.  Define 


P20  A*n 
Mu  M  02 


where 


L  =  y/pfiaU  S 
1  =  US2UT 


s/fct(Tj 

is  the  normalized  eigendecomposition  ofl.  Then  the  normalized  Frobenius-norm  error  is  given  by 


|Moo  kjsiLLT  —  X\\ p 

Wf 


KvMoo 


V^'nT  00 


Proof:  Letting  A  =  Moo kNLLT  and  B  =  X,  we  can  write 

A  =  aUS2UT, 
B  =  bUS2UT, 

where 


S2  = 


A  0 
0  1/A 

a  =  ki\rp  go, 


b  = 


^  AT  Moo 


(88) 

(89) 

(90) 

(91) 


(92) 

(93) 


(94) 

(95) 

(96) 
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epsilon 


Figure  21:  Relative  error  in  matching  second  order  moments  using  the  Initial  Guess  Algorithm 


Hence  we  have 


\\A-B\\f 


\\(a-b)S2\\F=\a 

\b\ 


b  I 


(97) 

(98) 


Hence, 

\\A- B\\F  _  \a-b\  _  VkNV oo  +  e  -  kNHoo  _  ^  fcjvMoo 

II-^IIf  H  VkNl4o  +  e  Vk'ifl4o  + e 

which  establishes  the  result. 


(99) 

□ 


We  have  plotted  the  expression  for  the  relative  error  in  Figure  21  for  N  =  3  and  N  =  1000,  and  assuming 
that  noo  =  1.  This  figure  shows  that  although  the  relative  error  grows  quite  fast  as  e  is  increased,  it  never 
exceeds  the  maximum  of  1  (i.e.  100  percent)  for  a  fixed  /zoo-  Also,  the  relative  errors  for  different  number  of 
sides  are  seen  to  be  very  close. 
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